### Reproduce TABLE 4 SRBI 

  rm(list = ls())
  setwd("")
  
  library(foreign)
  smoking <- read.dta("srbixx.dta")
  library(gam)
  library(Hmisc)

### Predicting Quitting: Relative vs. Attributable 
    smoking.gam1 <- gam(currsmok2 ~ s(rd) + s(rr),
                        subset = control==1&rd>-500,
                        data=smoking, family=gaussian(), weights=weight1)
    smoking.gam1a <- gam(currsmok2 ~ s(rd),
                         subset = control==1&rd>-500,
                         data=smoking, family=gaussian(), weights=weight1)
    
    smoking.gam1b <- gam(currsmok2 ~ s(rr),
                         subset = control==1&rd>-500,
                         data=smoking, family=gaussian(),weights=weight1)
  
  # LR Test from Adding Relative Risk to Attributable Risk: .03
    anova(smoking.gam1a, smoking.gam1)
    
  # LR Test from Adding Attributable Risk to Relative Risk: .64
    anova(smoking.gam1b, smoking.gam1)
  
  # Relative Coef: 13.8
    smoking2 <- smoking
    smoking2$rr <- wtd.quantile(smoking2$rr, weights=smoking$weight1, .25, na.rm=TRUE)
    smoking3 <- smoking
    smoking3$rr <- wtd.quantile(smoking3$rr, weights=smoking$weight1, .75, na.rm=TRUE)
    mean(predict(smoking.gam1, newdata=smoking2, type="response") - predict(smoking.gam1, newdata=smoking3, type="response"), na.rm=T)
    
  # Attributable Coef: 0.3
    smoking4 <- smoking
    smoking4$rd <- quantile(smoking$rd, .25, na.rm=T)
    smoking5 <- smoking
    smoking5$rd <- quantile(smoking$rd, .75, na.rm=T)
    mean(predict(smoking.gam1, newdata=smoking4, type="response") - predict(smoking.gam1, newdata=smoking5, type="response"), na.rm=T)

### Predicting Quitting: Relative vs. Absolute 

    smoking.gam1 <- gam(currsmok2 ~ s(pq52) + s(rr),
                        subset = control==1&rd>-500,
                        data=smoking, family=gaussian(),weights=weight1)
    
    smoking.gam1a <- gam(currsmok2 ~ s(pq52),
                         subset = control==1&rd>-500,
                         data=smoking, family=gaussian(),weights=weight1)
    
    smoking.gam1b <- gam(currsmok2 ~ s(rr),
                         subset = control==1&rd>-500,
                         data=smoking, family=gaussian(),weights=weight1)
    
    # LR Test from Adding Relative Risk to Absolute Risk: .002
    anova(smoking.gam1a, smoking.gam1)
    
    # LR Test from Adding Absolute Risk to Relative Risk: .15
    anova(smoking.gam1b, smoking.gam1)
    
    # Relative Coef: 15.2
    smoking2 <- smoking
    smoking2$rr <- wtd.quantile(smoking2$rr, weights=smoking$weight1, .25, na.rm=TRUE)
    smoking3 <- smoking
    smoking3$rr <- wtd.quantile(smoking3$rr, weights=smoking$weight1, .75, na.rm=TRUE)
    mean(predict(smoking.gam1, newdata=smoking2, type="response") - predict(smoking.gam1, newdata=smoking3, type="response"), na.rm=T)

    # Absolute Coef: -10.5
    smoking4 <- smoking
    smoking4$pq52 <- quantile(smoking$pq52, .25, na.rm=T)
    smoking5 <- smoking
    smoking5$pq52 <- quantile(smoking$pq52, .75, na.rm=T)
    mean(predict(smoking.gam1, newdata=smoking4, type="response") - predict(smoking.gam1, newdata=smoking5, type="response"), na.rm=T)

### Predicting Desire to Quit: Relative vs. Attributable 
    
    smoking.gam1 <- gam(dquit1 ~ s(rd) + s(rr),
                        subset = control==1&rd>-500,
                        data=smoking, family=gaussian(), weights=weight1)
    smoking.gam1a <- gam(dquit1 ~ s(rd),
                         subset = control==1&rd>-500,
                         data=smoking, family=gaussian(), weights=weight1)
    
    smoking.gam1b <- gam(dquit1 ~ s(rr),
                         subset = control==1&rd>-500,
                         data=smoking, family=gaussian(),weights=weight1)
    
    # LR Test from Adding Relative Risk to Attributable Risk: .09
    anova(smoking.gam1a, smoking.gam1)
    
    # LR Test from Adding Attributable Risk to Relative Risk: .27
    anova(smoking.gam1b, smoking.gam1)
    
    # Relative Coef: 17.0
    smoking2 <- smoking
    smoking2$rr <- wtd.quantile(smoking2$rr, weights=smoking$weight1, .25, na.rm=TRUE)
    smoking3 <- smoking
    smoking3$rr <- wtd.quantile(smoking3$rr, weights=smoking$weight1, .75, na.rm=TRUE)
    mean(predict(smoking.gam1, newdata=smoking3, type="response") - predict(smoking.gam1, newdata=smoking2, type="response"), na.rm=T)
    
    # Attributable Coef: -1.1
    smoking4 <- smoking
    smoking4$rd <- wtd.quantile(smoking4$rd, weights=smoking$weight1, .25, na.rm=TRUE)
    smoking5 <- smoking
    smoking5$rd <- wtd.quantile(smoking5$rd, weights=smoking$weight1, .75, na.rm=TRUE)
    mean(predict(smoking.gam1, newdata=smoking5, type="response") - predict(smoking.gam1, newdata=smoking4, type="response"), na.rm=T)
    
### Predicting Desire to Quit: Relative vs. Absolute 
    smoking.gam1 <- gam(dquit1 ~ s(pq52) + s(rr),
                        subset = control==1&rd>-500,
                        data=smoking, family=gaussian(),weights=weight1)
    
    smoking.gam1a <- gam(dquit1 ~ s(pq52),
                         subset = control==1&rd>-500,
                         data=smoking, family=gaussian(),weights=weight1)
    
    smoking.gam1b <- gam(dquit1 ~ s(rr),
                         subset = control==1&rd>-500,
                         data=smoking, family=gaussian(),weights=weight1)
    
    # LR Test from Adding Relative Risk to Absolute Risk: .05
    anova(smoking.gam1a, smoking.gam1)
    
    # LR Test from Adding Absolute Risk to Relative Risk: .06
    anova(smoking.gam1b, smoking.gam1)
    
    # Relative Coef: 13.9
    smoking2 <- smoking
    smoking2$rr <- wtd.quantile(smoking2$rr, weights=smoking$weight1, .25, na.rm=TRUE)
    smoking3 <- smoking
    smoking3$rr <- wtd.quantile(smoking3$rr, weights=smoking$weight1, .75, na.rm=TRUE)
    mean(predict(smoking.gam1, newdata=smoking3, type="response") - predict(smoking.gam1, newdata=smoking2, type="response"), na.rm=T)
    
    # Absolute Coef: -15.6
    smoking4 <- smoking
    smoking4$pq52 <- quantile(smoking$pq52, .25, na.rm=T)
    smoking5 <- smoking
    smoking5$pq52 <- quantile(smoking$pq52, .75, na.rm=T)
    mean(predict(smoking.gam1, newdata=smoking5, type="response") - predict(smoking.gam1, newdata=smoking4, type="response"), na.rm=T)
    